home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / KSINV.C < prev    next >
C/C++ Source or Header  |  1995-11-17  |  2KB  |  108 lines
  1. /* ============ */
  2. /* ksinv.c    */
  3. /* ============ */
  4. #include <math.h>
  5. #include <mconf.h>
  6. #include <miscdefs.h>
  7. #include <stdio.h>
  8. #include <stdlib.h>
  9. /* ==================================================== */
  10. /* LinTrp - Interpolates X0 Linearly in [X1,Y1],[X2,Y2] */
  11. /* ==================================================== */
  12. # if defined(__STDC__) || defined(__PROTO__)
  13. static double
  14. LinTrp(double X0, double X1, double Y1, double X2, double Y2)
  15. # else
  16. static double
  17. LinTrp(X0, X1, Y1, X2, Y2)
  18. double    X0, X1, Y1, X2, Y2;
  19. # endif
  20. {
  21.     return(X1 == X2) ? Y1 : Y1 + (Y2 - Y1) / (X2 - X1) * (X0 - X1);
  22. }
  23. /* ==================================================================== */
  24. /* KSmirInv - Functional inverse of KSmirnov distribution function    */
  25. /* ==================================================================== */
  26. # if defined(__STDC__) || defined(__PROTO__)
  27. double
  28. KSmirInv(int n, double p)
  29. # else
  30. double
  31. KSmirInv(n, p)
  32. int    n;
  33. double    p;
  34. # endif
  35. {
  36.     double  e, RetVal;
  37.     double  e0, e1, e2, p0, p1, p2;
  38.     /* ------------------------------------- */
  39.     /* Finds e such that KSmirnov(n, e) = p. */
  40.     /* ------------------------------------- */
  41.     if (p <= 0.0 || p >= 1.0)
  42.     {
  43.     if (p == 0.0)
  44.     {
  45.         RetVal = 0;
  46.     }
  47.     else if (p == 1.0)
  48.     {
  49.         RetVal = 1.0;
  50.     }
  51.     else
  52.     {
  53.         char   *ErrLabel;
  54.  
  55.         ErrLabel =
  56.         (p < 0) ? "KSmirInv(): p < 0" : "KSmirInv(): p > 1";
  57.  
  58.         mtherr(ErrLabel, DOMAIN);
  59.  
  60.         RetVal = -1;
  61.     }
  62.     }
  63.     else
  64.     {
  65.     /* ------------------------------------------------ */
  66.     /* Start with approximation p2 = 1 - exp(-2 n e^2). */
  67.     /* ------------------------------------------------ */
  68.     /* -------------------------------- */
  69.     /* Solve for e2 and Get Estimate p2 */
  70.     /* -------------------------------- */
  71.     e2 = sqrt(-log(1 - p) / (2.0 * n));
  72.     e2 = __min(e2, 1.0);
  73.     p2 = KSmirnov(n, e2);
  74.  
  75.     /* ------------------------ */
  76.     /* Get Estimate of (e1, p1) */
  77.     /* ------------------------ */
  78.     e1 = .9 * e2;
  79.     e1 = LinTrp(p, KSmirnov(n, e1), e1, p2, e2);
  80.     p1 = KSmirnov(n, e1);
  81.  
  82.     e0 = LinTrp(p, p1, e1, p2, e2);
  83.  
  84.     do
  85.     {
  86.         e = e0;
  87.  
  88.         p0 = KSmirnov(n, e0);
  89.  
  90.         if (p0 < p)
  91.         {
  92.         p1 = p0, e1 = e0;
  93.         }
  94.         else
  95.         {
  96.         p2 = p0, e2 = e0;
  97.         }
  98.  
  99.         e0 = LinTrp(p, p1, e1, p2, e2);
  100.     }
  101.     while (fabs(1 - e / e0) > 1e-10);
  102.  
  103.     RetVal = e0;
  104.     }
  105.  
  106.     return RetVal;
  107. }
  108.